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Abstract. We present an analytical study of a toy model for shear banding, without normal stresses, which 
uses a piecewise linear approximation to the flow curve (shear stress as a function of shear rate) . This model 
exhibits multiple stationary states, one of which is linearly stable against general two-dimensional pertur- 
bations. This is in contrast to analogous results for the Johnson-Segalman model, which includes normal 
stresses, and which has been reported to be linearly unstable for general two-dimensional perturbations. 
This strongly suggests that the linear instabilities found in the Johnson-Segalman can be attributed to 
normal stress effects. 

PACS. 47.50.-d Non-Newtonian fluid flows - 47.20.-k Flow instabilities - 47.57.Ng Polymers and polymer 
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1 Introduction 

Shear-banding phenomena have received a lot of attention during the last decade. In shear banding the material splits 
into different spatial regions, or bands, that flow at different shear rates 7. This phenomenon has been observed in 
granular media [1] and in viscoelastic living polymer systems such as wormlike micelle solutions [2]. Theoretically, 
shear banding is fairly well understood at the level of a stable one dimensional (ID) banding profile [3l|4l[5j, which 
connects two shear rates that are stable at a given selected total shear stress E. Wormlike micelles, polymer melts, and 
liquid crystals naturally give rise to such bistable constitutive relations [6l[7l[8]. However, the stability of bands in two 
or three dimensions (i.e. with respect to capillary-like fluctuations of the interface) is not well understood, despite a 
frequent interpretation in terms of a simple stable flat interface between coexisting states. It has been experimentally 
observed that the interface between the two phases is not necessarily flat, but can exhibit strong undulations and 
erratic fluctuations [9l[T0llllJ. Theoretical calculations have shown that in shear banding flows, chaotic motion of shear 
bands can occur jl^lllSliTl] . Such chaotic motion has also been inferred experimentally in recent works by Sood and 
co-workers [T5l[T6|[T7] . 

A recent numerical calculation in two dimensions, which was first performed by Fielding [18] and whose results were 
later confirmed in Refs. 19, 20 , demonstrated that, for the Johnson-Segalman model the interface between the low 
and high shear rate phases is linearly unstable to undulations. In this case the unstable mode involved normal stresses; 
however, it is still not known whether normal stresses are carried by an instability inherent in the two-dimensional 
(2D) nature of the fluctuation, or whether normal stresses trigger the instability |20| . Hence, in this paper we consider 
a simple general toy model without normal stresses [4lll4|. We show that for a wide class of multivalued flow curves, a 
linear instability can never occur. Our findings imply that the coupling between convective terms and perturbations 
in shear stress cannot lead to a linear instability, which suggests that an instability requires other degrees of freedom, 
such as normal stresses. 

This paper is organised as follows. In Section [2] we introduce the model and study stationary solutions for general 
constitutive equations in which the shear rate is multivalued for a range of stresses. 2D perturbations are introduced 
in the model in Section [3] and a linear stability analysis is carried out, from which we conclude that linear instability 
docs not occur for our system. This strongly suggests that normal stresses are responsible for the linear instability 
observed in Refs. [18]. The results are discussed and summarized in Section HI 



2 Model Description 

We consider planar Couette shear flow between flat plates separated by a distance /i, for which the velocity field is 
V = jyx. In the very low Reynolds number limit, which applies to the complex fluids of interest, the total shear stress 
S is uniform in space. We assume that S comprises two terms, 

S^ap{y) + l^jiy), (2.1) 

where the second term is the Newtonian stress of the solvent and a-p is the polymer stress. In this paper, we only 
consider the shear component of the stress tensor, but take the 2D nature of the perturbations in the fiow field into 
account. Following |14(f3P4]. we consider the following governing equation for the polymer shear stress: 

(9t + V V) = + -.g(7r) + D^^Up, (2.2) 

T T 

where r is the relaxation time of the polymer stress, D is the stress diffusion coefficient, and G is the plateau modulus. 
The function g is nonmonotonic, and in Refs. [¥1111] was taken to be 

Stress diffusion is necessary to define a uniquely selected stress S* [5], as is usually found in experiments under 
controlled shear rate conditions. Moreover, it is natural that spatial gradients in the microstructurc should be penalized, 
no matter how weakly. 

We next make all quantities dimcnsionlcss, as noted by a carat, by expressing stress relative to G, length relative 
to the plate separation h, and time relative to the relaxation time r: 

ap = a/G S = E/G (2.4) 

= r7 D = Dr/h? (2.5) 

y — yjh X = x/h (2-6) 

i=t/T v^vr/h (2.7) 
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Fig. 1. The polymer stress (jp as a function of the distance y = y/h from the lower plate for D = 0.0001 (dotted), D — 0.001 
(dashed) and D = 0.01 (solid), for a — 20 and average shear rate {7) — 2.60. The selected stress is E* = 0.470 for D = 0.001 
and D — 0.0001; and S* = 0.471 for D = 0.01. The inset shows the corresponding constitutive curve, shear stress X" as a 
function of 7, for a = 20. 



The total stress is a nonmonotonic function of shear rate 7. For an average shear rate (7) imposed on the unstable 
portion of the flow curve, with negative slope dS /dj < 0, the system will break up into two bands, and hence a spatial 
variation along y. Upon integrating Eq. (|2.ip over the vertical distance between the plates we find 



where spatial averages are defined by 



a 



{O) = / 0{y)dy 
Jo 



(2. 



(2.9) 



and we have introduced the dimensionless parameter where a — Gt/tj. For an inhomogeneous steady state, the polymer 
stress equation, Eq. (|2.2p . leads to 

- ap + gia{S - ap)) + m; ^ 0, (2.10) 

where the strain rate has been eliminated using Eq. (j2.ip and the prime denotes a y-derivative, a'p = dup/ dy. As in 
[T8l[2T|. we impose Neumann boundary conditions, (y'p{y = 0) = (y'p{y = 1) = 0, for which Eq. (j2.10p can be solved 
numerically for a given stress S. 

The selected stress S* is determined by the constraint of fixed average shear rate (7) . We thus eliminate S from 
the stationary condition, Eq. (|2.10p . using Eq. (|2.8|) leading to 

- CTp + 5 [a ({Gp) - dp) + ('^)] + ba'l, = 0. (2.11) 

We numerically solve this integro-differential equation using a semi-implicit Crank-Nicolson algorithm [22j. Once a 
solution is found, the selected stress is obtained using Eq. (|2.8p . Fig.[T]shows banding profiles (Tp{jj) for the function 
g{^) given by Eq. ((23)) . 

To better understand stress selection we will generalize g{S^) to a function that allows explicit analytic calculations, 
but is sufficiently general to describe shear banding. The calculations discussed in the main text concern this general 
model. However, for comparison we also reproduce the calculations for the function g defined in Eq. (|2.3[) . in Appendix 
ICl The most essential characteristic of g is that it vanishes for small and large shear rates, so that the function 
S{x) = x/a + g{x) is nonmonotonic with a maximum and a minimum. A simple piecewise linear function showing 
this behavior is 



^0(7) 



a 

Acio - 71 



70 - 71 



(7-70) ^c7o 



7 
V. a 



(7 < 70) 



(70 < 7 < 71) 



(7 > 71), 



(2.12) 
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Fig. 2. The selected stress S* and the piecewise linear function X'o(7) as a function of 7 for 70 = 0.2, 71 — 0.8, Ac — 5; the 
selected stress is S* — 0.891. 
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Fig. 3. The solutions A, B, C, for yo and yi, for the parameter values: Ac — 5,a = 1, (7) — 0.50, 70 ~ 0.2, 71 = 0.8, and 
D = 0.001. The dashed curves correspond to Eq. (|2.17|) and the solid ones to Eq. (|2.18fl . The physical range of solutions is in 
the upper left triangle, < < J/i < 1- 

where Ac'jo > 71. The form of ^0(7) is completely specified by the three parameters 70, 71, and Ac- 
Using -£'0(7) ~ ilea instead of g in the equation of motion for (jp, and eliminating using Eq. (|2.1[) . the steady 
state condition for 7 is given by 

hi" + ai: ^ al^^ii). (2.13) 
This is solved vifith Neumann boundary conditions for the stress at y = and 2; = 1, leading to 



7/(y) = ci cosh { j + 



in{y) = cos ( \l^y) + sin (j 



iO<y< yo) 



- 70 — (70 - 71) 

^c7o - 71 



iiii{y) = OiS + C2 cosh 



y-1 



{vi<y< 1) 



(2.14a) 



{yo<y< yi) (2.i4b) 



(2.14c) 
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Fig. 4. Two stationary solutions corresponding to A (solid curve) and B (dashed curve). An additional third solution with one 
more oscillation is not shown. 



where z = {A^'jo — 7i)/(7i — 70)1 which can be recognized as the negative of the slope of the function £"0(7) between 
70 and 71. The profile obeys the following continuity conditions: 

iiivo) = jiiim) = 70 (2.15) 

inim) = iiiiiiii) 71 (2-16) 
^iiyo) = ^ii{yo) (2.17) 
iiAyi)=^iuiyil (2.18) 



where again the prime denotes a y-derivative, 7' ~ dj/dy. These six conditions guarantee a smooth solution for the 
shear rate and stress and from them the four constants ci, C2, di, ^2 and the locations yo and yi where the function 
So{j) switches from one piece to the next, can be determined. Appendix [Al contains expressions for the constants ci, 
C2, and di, d2 as functions of yo and yi. di and d2 are found by combining Eqs. (|2.15p and (|2.17p . Eqs. (|2.16p and 
(|2.18p are then used to determine the values of yo and yi . The selected stress S* is determined by the constraint of 
an apphed average shear rate (7) , which is given by Eq. (|A.2[) . In a shear banding state the selected stress S* satisfies 
71 < 17* < Ac7o, and we must have 7 > on the interval y G [0, 1]. 

A graphical solution for yo and yi, for the case D — 0.001, Ac ~ h and (7) = 0.5, is depicted in Fig. [3l as the 
intersections of the contour plots of Eqs. (|2.17p (dashed hne) and (|2.18p (solid line), which should simultaneously be 
obeyed by yo and yi. To solve these equations the expressions for the constants ci, €2,^1,^2 from Appendix 1X1 are 
used, Eqs. (|A.1[) . as well as expression (|A.2[) for the selected stress. The physical range of Fig. [3] is the upper left hand 
triangle < yo < yi < 1, to the left of the diagonal line in Fig. [3l 

Fig. [3] shows three intersection points corresponding to stationary solutions, denoted hy A = (0.4806,0.6101), 
B = (0.2984,0.7726), and C = (0.1160,0.9345). The values of the selected stress S* associated with these solutions 
are 0.8944 for A, 0.8946 for B and 0.8949 for C . The fourth intersection point (the companion oi A, with the unphysical 
value of S* = 4.8) is ignored, as it violates the condition 71 < < Ac'^q. 

The first solution A corresponds to the common shear rate profile depicted as a solid curve in Fig. |4l The stationary 
solution B (shown dashed) has an extra oscillation in the interface region (Fig. [4|) , while the third solution C has two 
oscillations (not shown) . The selected stress changes monotonically as a function of an increasing number of oscillations 
in the resulting solution. The existence of multiple solutions for small diffusion constants is well-known in the context 
of pattern formation in reaction-diffusion models [53]. Multiple solutions occur in our model when D < 10^^, with the 
multiplicity increasing for decreasing D. Having obtained the stationary solutions we next examine their stability by 
performing a linear perturbation analysis. 
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3 Linear stability analysis 

3.1 Governing Equations and Matching Conditions 

We study the linear stability of the piecewise linear model by introducing small perturbations to the stationary ID 
solutions, 

0'p(i,y) = CTp^(y) +(S(Tp(x,y) (3.1) 

S{x) ^ + 6S{x) (3.2) 
'j{x,y) ^ j^^{y) + Sj{x,y) (3.3) 
V = (v^^iy) + Sviix,y),Svy{x,y)), (3.4) 

where the superscript denotes the stationary ID solution. Perturbations in the shear rate, polymer stress and total 
stress are related by SU = Sup + 7/0;, which can be recast as 

{5ap) ^ 5&P + 5-^ /a (3.5) 

by averaging and using ((57) = 0. Note that {5'^) ~ implies that the average perturbation in the polymer stress 
{6ap) equals the perturbation in the total stress SS, for 7 — a(6S — Sap). The shear rate is given by 7(2;, y) = 
dy{v^^{y) + dvx{x,y)), so that 5'j{x,y) = dy5vx{x,y). Because the flow is incompressible, a variation of the velocity 
in the x-direction must be accompanied by a compensating variation of the velocity in the y-direction. which obeys 

dy5vy{x,y) ^ ~di5vx{x,y). (3.6) 

We recall that the equation of motion of the polymer stress is given by 

{di + w-V)ap = SQ{^l)~S + DV^-CTp. (3.7) 

This equation can be recast in a differential equation for the shear rate 7, by using Sup = 5S — 5'^ /ol^ from which we 
find to first order in all perturbations 



+ b [dxi + dyy] - bdxiaSS - aSoiV + S^) + aSoiV), (3.8) 

where we explicitly kept £'0(7^^ + ^7) and So{j^^) in order to avoid differentiating the function Sq with respect to 
7. This was done as the function So is continuous for all 7, but non-differentiable in the turning points of the stress 
at7o,7i. 

Using Eqs. (|3.5p and (|3.6p we can easily express the perturbed velocities Svx{x,y) and Svy(x,y) in terms of Sj. 
Starting from the identity ^^^^^■y^^'^ ~ ,5^^ ^^jq integrate from to y and use the no slip boundary condition at y = 0, 
which yields 



5vx{x,y) 



5i{i,y')dy'. 



(3.9) 



If one next differentiates Svx{x,y) with respect to x and substitutes the resulting expression in Eq. (|3.6p . only one 
integration with respect to y (again starting from and using no slip boundary conditions) is necessary to obtain the 
perturbed velocity in the y direction 



Svy{x,y) = - 



dxSj{x,y")dy" 



dy'. 



(3.10) 



Hence the equation for the shear rate perturbation Sj (|3.8p can be expressed entirely in terms of Sj, integrals over Sj 

and SS. Perturbations in the shear rate must be continuous across the matching points of the bands, which are yet 
to be found. These follow from demanding 



yo+Sya 



7o> 



^7 + (57 



vi+Syi 



7i> 



(3.11) 
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that is, the matching points arc shifted by perturbations 6y needed to bring the local shear rate to the matching points 
in the constitutive relation 170(7). Hence we split the interval [0, 1] into three parts: (/), y g [0,ya]: {II), y € [yo,yi], 
and {III), y £ [j/i, 1], where jji = yi + 5yi{x, t), i = 1,2. The equations of motion for the perturbations in the polymer 
stress are, from Eq. 



di5^^' = DW^6Y - vl''{y)ds,5y - Svy{x, y)dy^^'' - AXi' + F{adS, t) (3.12a) 

diS'Y' = D^^Sl" - vl^{y)d±6Y' - 5vy{x, y)dy^^^ + zS^' + F{a6S, t) (3.12b) 

d-Xi''' = DV^S^i''' - vl°{y)di5'^'" - 5vy{x, mii^"" - S^''' + F{aSS, t) (3.12c) 

where we defined F{a5E,t) = {d^ + v\^dx — Ddxx)a5S{x,i). Eqs. p.l2p must be supplemented with boundary 
conditions the matching points. At y ~ i/Q this condition can be found from Eq. (|3.1ip . 

¥^{yo + Syo{£,i)) + S^yo + 6yo{x,i)) ^ ^0, (3.13) 
which leads, to first order in the perturbed quantities, to 

SHyo,x,i) = ~l^^] dyo{x,i). (3.14) 



dy 



Vo 



Eq. (|3.14p relates the change in the shear rate to a shift in the interface position. The boundary condition at y = yi 
is derived similarly. From Eq. p.l4p we immediately retrieve the shift in the value of yo, and hence the new position 
of the two shear bands, once S'^ is known. To deal with the inhomogeneous (time dependent) boundary conditions in 
region (//) we define a new function q{x,y,t) by 

,(i, y, i) = ^^-(i, y, i) + iy-y^)^^^yo^-^i) _ {V'V.mm,S:,i) ^ „^ ^ 

yi -yo 2/1 - yo 

which satisfies continuity at the matching points. Hence, the function q{x,y,t) vanishes at the boundary points by 
construction. In Eq. p.lSp we introduced the function G{x,y,t) for notational convenience. Using Eq. (|3.12bp . we find 
the following differential equation for q: 

diq ^zq + L>V\ ~zG + F{a5S,i) + F{G,i) - vl"{y)d£q - 6vy{x,y)dy^^° . (3.16) 

Now that we have derived the governing equations for the perturbations, we study their stability properties by 
employing a Fourier expansion in spatial coordinates. 

3.2 Fourier Expansion 

We next expand the shear rate perturbation S'j in Fourier modes, within the three regions, consistent with the boundary 
conditions 7' = at j) = and y = I. This results in a continuous but non-smooth function 6j over the entire interval 
[0,1]: 

dY{x, y, i) = V cos (^^^i^i±ll) A„(£, t) + <5^(yo, x, t) (0 < y < yo) (3.17a) 

;s V 2yo ; 

g(i,y,i) = Vi3„(x,t)sinf^^^fc^') (yo < y < yi) (3.17b) 

5^^''\x, y, i) ^ { '^^^2(1^^'^) ) + ^' ^ (yi<y<l)- (3-17C) 

Only odd cosine modes contribute in regions / and /// because we require zero derivatives at the ends of the interval 
[0, 1] and a matching condition at the right, respectively left, ends of the intervals [0,yo] a-^d [yi, 1]. If one substitutes 
y = yo in Eq. (j3.17ap one finds 5'^^ = '57(yoj x,t), which coincides with the boundary value of 5'y in yo, derived above 
in Eq. (|3.1ip . which guarantees contimuty of S'y{x,y,t) throughout the interval [0, 1] 
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One can easily demonstrate that the solution of Eqs. ()3.17p is non-smooth, by trying to match the derivatives 
at y = yo, yi using Eqs. (|3.17p . We will not reproduce these calculation here, but only note that the narrower the 
interface, that is the smaller yi — y^, the closer the matching at the the interface points j/o and j/i. The piecewise linear 
model will in this limit approach the smooth model defined by Eq. (j2.3p . 

We therefore conclude that the perturbations are in general only continuous. To obtain an equation linking 
Sj{yo,x,t) and S'j{yi,x,t) we invoke the constraint {6j) = 0, from which we find 



[Vo + yi) + 



(S7(yi,x,£) 



(2 - 2/1 - yo) 



E 



A„(-l)"2yo 



E 



2(l-yi)A.(-l)^ 



(2n+l)7r ^ (2n-fl)7r 

n=0 ^ ' ri=0 ^ ' 



(2/1 -yo)E 



(3.18) 



Eq. p.lSp constitutes the final equation needed to close the system, in a way which is similar to the method we used 
in section [2] to deduce the value oi as we will now demonstrate. 

The general idea is to find the governing equation for the amplitudes A„, i3„, £)„ from Eqs. (j3.12p by projecting 
out the cosine (for An and D„) and sine modes (for -B„). The evolution for 5'^{yi^ x, t) and Sj^ya, x, t) can be obtained 
by integrating Eqs. (|3.12ap and (j3.12cp over their respective intervals. This is most easily done in the Fourier domain, 
so we expand An{x,i), Bn{x,t), Dn{x,t), SIJ{x,t), and S'j{yo,x,t) and 5j{yi,x,i) in Fourier transforms: 



An{x,i)= an{k,i)e'^^ dk, B^ixJ) = bn{k,i)e'''^ dk, D^ixJ) = dn{k,i)e'''^ dk, 

J — oo <J —oo J —oo 

nOO nOO nOO 

SS{x,i)= 5S{k,t)e'^^ dk, S^{yo,xJ)^ S^{yo,k,i)e'''^ dk, S^{yi,xJ)= S^{yi,k,i)e'''^ dk. {3A9) 



From Eqs. (|3.12ap and ()3.17ap we then obtain the following evolution equation for the amplitudes an{k,i) and 
6'y{yo,k,t) The amplitudes an{k,i) can be extracted from this equation by projecting out the different cosine modes, 
whereas direct integration gives the evolution equation for 5j{yo, k,t). We can perform the same calculations on 6j^^^ , 
and q which leads to similar expressions involving dn{k,i), 6'y{yi, k,t), and bn{k, i), which are derived in Appendix IB] 
It turns out that the introduction of scaled amplitudes («„, /?„, i5„) is beneficial for numerical and notational reasons; 
they and related to (a„, 6„, dn) according to 



a„(-l)"(2n+ l)7r 



hr,n-K 



<i„(-l)"(2n-h l)7r 



(3.20) 



To show the form of the evolution equation of the amplitudes, we reproduce the governing equation for the scaled 
amplitudes Q;„(A:,t) 



^7 



ID 



dam 
di \ dy 

+ bk^aS{k,i) 

OO 

+ E 



di5yo{k, t) - dia5S(k, t) = {A^ + Dk') 



Vo 



dy 



Sya{k,i)- {Ac + Dk'' + D 



yo 



7r(2?7i+ 1) 



2yo 



ifc(-l)'"(2m-t-l)7r 



E""(-i)'^ 



16yo rpA ^Ai , / di 



ID 



n=l 



4a„(-l)" (jA^D_j.A 



(2n+ l)37r3 L ™" 



dyo 



5yo{k.i)F^ 



{2n+ l)7r 



Ar 



dy 



5yoik,t) 



(3.21) 



The quantities F^f^, R^m defined in terms of integrals over the interval [0, yo\. Their precise definitions 

can be found in Appendix IB| Eqs. (|B.3p . 



Finally, we need one further differential equation for the evolution of 5E{x,i). This equation is found by dif- 
ferentiating the identity (|3.18p with respect to time and subsequently substituting the expressions for dfd'jitjo, k,t) 
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. 4 



yo 



Fig. 5. Comparison of the analytical solution for yo and yi, depicted by the intersection of the two solid curves given by 
(Eq. (|3.27p ) and (Eq. (I3.24p ). with the solution obtained by the method of section[2l which are printed dashed. The parameter 
values are Ac = 5.0, 70 = 0.2, 71 = 0.8, {7) = 0.5 and D = 0.01. The diagonal curve indicates the line yi = yo and physical 
solutions are in the upper left triangle. The intersection point A between the two dashed curves is practically on top of the 
intersection point of the solid curves {yo = 0.3033,yi = 0.7130). 



(Eq. (|R2l)) and di5^{yi,k,t) (Eq. (|R5)) . This yields 



00 J - 00 jr. 



^ dt (2n-hl)27r2 

00 J- 

+ 4(l + I)fc^)(2-^,-^o)E(^+l)... 



dt (2n-hl)27r2 

n—0 ^ ^ 



H-4(A, + m2)(yi+y„)E 



of ' \ d(32n+l 

8(yi-?/o)2^ 



n=0 



di {2n + l)27r2 



Tl = 



n=0 ' 



D{2-yo~yi] 



„ . -J^^ + (l + m-2)(2-yi-yo)^ + (A, + m2)(yi-fyo)^-«i3fc25£ + ifc77(fc,t), 

(3.22) 



where n{k,t) contains the imaginary contribution to (j3.22p . whose exact form can be found in Appendix [Bl 

To summarize, the linear stabihty of the shear banding states to undulating perturbation is investigated by an- 
alyzing the evolution of the amplitudes Q;„(fc,i), f3n{k,t), (5„(fc,t), SIJ{k,t) and S'y{yo, x,t), Sj{yi, x,t) governed by 
Eqs. (|3?2T|) . (jB.ep . (|B.4jl . (|B.2p . (jB.SP and ([3221) • The profile, and hence the interface, is stable only if all these 
amplitudes decay with time. 



3.3 Non-undulating perturbations for D^O 



Before turning to the general case, we will examine the limit of vanishing D and k ~ 0. This special case was 
investigated for the Johnson-Segalman model in [13] and was also explored in The conclusion in both cases was 
that the shear bands in the system are neutrally stable in the absence of diffusion terms. 

In section [5] we demonstrated how to find the stationary state(s) of system defined by Eq. (|2.2p graphically as the 
intersection point (s) of two families of curves, defined by Eqs. (|2.17p and (|2.18p . For the case D ^ 0, we can actually 
find an excellent analytic approximation to the stationary solution(s) which is remarkably close to stationary solutions 
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50 100 150 200 

smallest 200 eigenvalues 

Fig. 6. Real parts of the 200 eigenvalues with smallest real parts of the system of equations (|3.28|) for non-undulating pertur- 
bations. The parameter values are Ac = 5,yo = 0.4806, yi = 0.6101, (7) = 0.5, 70 = 0.2, 71 = 0.8, and D = 10~^. The smallest 
eigenvalue is 0, demonstrating that the system is neutrally stable. Eigenvalues are more concentrated around A = 1 and A = 5 
in agreement with the analytical asymptotic result. The first eigenvalue A = is related to SS, the second with A ~ 0.1 to /3i 
and all eigenvalues A = 1 to Qn , and \ = Ac = 5 to Sn- 



found by the graphical method. In the hmit D ^ equations (|2.16p and (|2.18p reduce to 
aS - Ac^o 



cos I w4(yi - yo)^ ^ ^ X") V (\/^*^^^ " ^^"O " \ ^ '~X^i~~^] ^ ^"^ (3.23a) 



aS - Acjo 



-{yi -yo)j = 7= — , 

D Jl 



(3.23b) 



keeping all terms up to order V D. If we next square Eqs. p.23ap and (j3.23bp and add them we find the following 
relation 



aZ'-Ac7o\ Ac(z + 1) 



ai7 — 71 



(3.24) 



If we would next substitute the value of S as given by Eq. (|A.2p . we could express y\ in terms of yo- However, first we 
obtain the width of the interface, that is the value of yi — yg, to verify if it scales proportionally to in agreement 
with findings in the literature. We therefore solve Eqs. p.23p for sin {^\J^{yi ~ yo)j ^ind obtain, after using Eq. p.24p . 



sni I , / -T-(yi - ijo) = , ^ ^ , ■ 

D J ^z+l+A^ + AJz 



(3.25) 



When we invert this equation we get the following expression for interface width 



yi-yo = \ -( (2A/ + l)7r - arcsin ( ^ = ) ) • (3-26) 

V^V \^/z + l + Ac + Ac/z J J 

The other solution of Eq. p.25p is not compatible with the condition 71 < aS < AcJq and therefore does not 
correspond to a physical solution. The index M is the band index. M — designates the two band solution (A), 

M = 1 has two extra bands, like solution B, and so on. Expression p.26p confirms that the interface scales as \/d, 
as observed in [UIIl], and gives an exact value for the interface length as a function of z and Ac- An implication of 
Eq. p.26p is that shear bands can only exist when D is sufficiently small. Indeed the maximum possible width of the 
interface is yi — j/o — Ij implying a maximum value of D* « (tt-i)'^ ■ case we have been considering z ~ 1/3, 
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— ' 1 

* solution A 

• solution B 
■ solution C 



Fig. 



smallest 10 eigenvalues 
7. The real part of the first 10 eigenvalues of S~^L'^ for an expansion of Eq. 



(|2.1|l in n = 1000 modes with zero wave 
number k for solutions A (A), B (o), and C (□). The parameters are the same as in Fig. 3: Ac = 5, (7) = 0.50, and 70 = 0.2, 
71 = 0.8; the diffusion constant D — 0.001. It can be clearly seen that the real part of the eigenvalues for solutions B and C 
cross zero, indicating instability of the solutions with multiple bands. Only the solution with two bands (A) is (neutrally) stable. 



and Ac = 5, we find D* ~ 0.054, which agrees with our numerical findings. The values of yo and yi for small D can be 
determined analytically by using Eq. (|3.24p once more, and substituting the expression of aE* as given in Eq. (jA.2[) 
thereby using Eq. (|3.26p . This calculation gives the following expression for yo, which designates the position of the 
interface of the steady state with 2(M + 1) shear bands 



(^c-l)(V^c7o + V71 
7o7i 



(71 - (7)) V Acio + (^c7o - (7)) V 71 



D 1 



(2M +l)TT-h 

7i 



(Ac 



1)^ 
(3.27) 



where h = arcsin ( VA^+i \ ^ q jj f^j. ^j^g parameter values used in this paper. From Fig. [5l one can see 

that the analytical solution (|3.27p with M = found for I? ^ 0, is still an excellent approximation for D = 0.01. 
The other stationary states with multiple {2{M + 1)) bands can be found similarly, by substituting the corresponding 
value of M in Eq. p.27|) . Note that multiple bands are not present for D = 0.01, as is shown in Fig. [5] where we only 
find one intersection point (A) of the solid curves. 



In order to study the stability of the solutions obeying Eq. p.26p . to non-imdulating perturbations we need to 
determine the eigenvalues of the evolution equations for a„{t), bn{t), dn{t) (or equivalently a„(t), /3n(t), (5„(t)), Sjijjo, t), 
^i{yo:t) and 6IJ{i) for fc = 0. Because the imaginary parts of the governing equations all depend linearly on fc, we 
can restrict ourselves to the real part of the differential equations for the amplitudes a„, /3„, (5„, Sj{yQ,t), S'j{yo,t) 
and SS. These equations can easily be found by projecting out the different sine and cosine modes and are given by 
Eqs. jBT]), ((R6l) . dEU, ((R2|) . (|B3|, and (|R9|) . For fc = and D the evolution equations for a„(fc,£), /3„(fc,£), 
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Sn{k,i), Sj{yoJ), Sj{yi,i) and 5S{k,i) read: 



dam 
di 

df3,n 
di 



dads d5'y(yQ,i) 



dt 



dt 



daSS /I - (-1) 



dt 



2 

z(-l)" 



-AcSj{yo,i) - Acam 
IdS^iijoJ) i-ir d6^{yi,i) 



dt 



d6m dads dSj{yi,t) 
di di 

dSj{yQ,t) daSS 

di di 

d6j{yi,t) daSS 

di di 

daSS 



dt 



-Siiyi,t) - 6„ 
dan 8 



di (2n-fl)27r2 

n— ^ ^ 



dt 



1- 



{{2M + i)TT-hy 



I3„ 



E 



-AJ^{ya,i) ~ A^Y^ 



E 



dSn 



di (2n-hl)27r2 

71 — ^ ' 



= -SiiyiJ) - ^ 



(2n+l)27r2 

8Sn 



(2n-h l)27r2 



n=0 



dt 



A.Xyo + m)'-^ + i2-y, 



yo) 



4A,(yo + yi)^ 



(2n + l)27r2 



+ 4(2 - yi - yo) 



ri=0 



(3.28a) 

(3.28b) 
(3.28c) 
(3.28d) 

(3.28e) 

Sn 



(2n + l)27r2' 
(3.28f) 



Eqs. p.28p contain a lot of information. It is immediately clear, for example, that A = is an eigenvalue of Eqs. p.28p . 
by which we mean that solutions with all amplitudes decaying as exp(— At) with A = 0, exist. This can be seen by 
interpreting Eqs. p.28p as a matrix equation which has the eigenvalue corresponding to 5S equal to zero. 

If it is proved that all amplitudes except = are always decay, and thus have eigenvalues with positive real 
parts , we have shown that the shear banded state is neutrally stable. This agrees with Fielding [T3] and Yih |24j . 
It should also be noticed that the equation for SS decouples from the other equations as no SS terms appear in the 
equations for the other amplitudes and we could therefore climate all dfSS terms in Eqs. (|3.28p by substituting the 
right-hand side of Eq. (|3.28f|l for dfSS in first five equations of (|3.28p . Thus we only need to focus on Eqs. (I3.28all3.28el) . 



By subtracting Eqs. (|3.28dp from p.28ap . we can find an equation expressed entirely in terms of a„ . and subtracting 
Eqs. (|3.28ep from p.28cp gives an equation in terms of (5„: 



dam 
di 



^E 

n=0 
oo 



dan 



1 



dSm g 



di (2n-Hl)27r2 

d^n 1 

"dT (2n l)27r2 ' 



~Ar_ 



a,, 



^E 



(2n-t- l)27r2 



Sn 



4 {2n + 1)2^2 

n—0 ^ ' 



(3.29a) 
(3.29b) 



The eigenvalues of Eq. p.29ap arc \ ^ Ac and of Eq. (|3.29bp A = 1, reflecting the decay of these amplitudes in time. 

When we next consider Eqs. (|B.2p and (|B.5P and Eq. p.28f1 to eliminate the terms, we find that the perturbations 

in the interface shear stresses obey 

di 

c% _ 
di 



Ac [l- 

yo + m 



m + yi 

2 

^i + A 



(^70 



yo + yi 



yo + yi 



70- 



(3.30a) 
(3.30b) 



Using Eq. (|3.18p with the amplitudes An and Dn set to zero, as we have previously found that these decay in time, and 
neglecting terms of order yi — yo we readily find that Eqs. (|3.30ap and (|3.30bp reduce to a single differential equation 
for 7o that reads 



dSjo 
di 



-Ac (1 - yo) Sjo - Vo&io, 



(3.31) 



which clearly corresponds to a solution that decreases in time. As in this approximation 5'yi is proportional to 570 the 
same holds true for ^71, showing that the the position of the interface is linearly stable. This means that the stability 
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of Eqs. (|3.28p only depends on the properties of the amplitudes Pm- The even and odd modes of the amplitudes Pm 
satisfy different evolution equations. Using the fact that a,i, (5„ and Sjo, ^ji decay in time, the stability properties of 
the amplitudes /?„ for both odd and even n are determined by the equation 

di ((2Af-f l)7r-/i)2 

where M is the band index. The system of equations (j3.32p gives eigenvalues 

As the smallest eigenvalue is Ai , the sign of this quantity determines the stability of the amplitude of /3„ and therefore 
the stability of the system. For M = 0, Ai > 0, which proves that the two banded solution is neutrally stable: all 
amplitudes a„, /?„, (S„, Sjo, Sji are decaying in time and a6E is staying constant in time. Multiple banded modes 
with M = 1, 2 all have an eigenvalue Ai < 0, which indicates instability of these modes. 

In conclusion we have proved that for small D the stationary shear band solution, with two bands is neutrally 
stable against non undulating perturbations and all other stationary states are unstable. Our analytical findings can 
be verified numerically. In Fig. [H] the results of a numerical calculation of the 200 eigenvalues with smallest real part 
are displayed, for D ~ 10^^. Besides the zero eigenvalue corresponding to SS, we indeed find X = Ac corresponding to 
the amplitudes a„, and A = 1 corresponding to the amplitudes 5„, many times. The degeneracy is lifted by the small, 
but finite value of D and the finite size of the matrix. The eigenvalue A « 0.1 corresponds to To check our result 
that all multiple bands solutions are unstable, we numerically calculated the smallest eigenvalues of Eqs. ()3.28|) for 
the solutions B and C, found in section[2l We indeed find that solutions B and C are unstable to small nonundulating 
perturbations. This is shown in the graph in Fig. [71 where an eigenvalue with real part smaller than zero was found 
for both solutions B and C. 

We next turn to the general case in which the diffusion is no longer assumed to vanish and undulations in the 
perturbations are admitted. 



(3.32) 



3.4 Linear stability: nonzero difFusion 

The linear stability of the stationary solutions found in section[2]with respect to undulatory perturbations are examined 
numerically. We follow the same strategy as in the previous section, projecting out the different sine and cosine modes, 
but this time keeping the imaginary parts of the evolution equations (|3.17p as k is no longer presumed to vanish. 

It again proves beneficial to use the scaled amplitudes Q!„(A:,i), /3„(fc,t), 5n(k,t) defined in Eq. (j3.20p rather than 
the amplitudes an{k,t), dn{k,t), for numerical convenience. Once we have projected out the sine and cosine 

modes, we are left with the evolution equations of the amplitudes. The details of this calculation are relegated to 
Appendix [Bl Here we merely note that the general structure of the system of differential equations, whose stabihty 
we would like to explore can easily be captured in a matrix equation. 

By introducing a vector u = {ao, ■ ■ ■ , q;„, /3i, • • • , /3„, Sq, ■ ■ ■ , S-jOi ^7i, S^), the matrix equivalent of the ampli- 
tude evolution is given by 

S ■ — = -CL^ - 2ikL^) ■ u, (3.34) 

dt 

where the (3n + 5) x (3n + 5) matrices S and L^,L' can be easily read off from the amplitude equations (|B.mB.4|) 
in Appendix IB] As in the previous section eigenvalues of the matrix S^^L = S^^L^ — 2ikS^^lJ- with positive real 
part correspond to stable solutions. The linear stability of the stationary solutions found in section [3.31 with respect 
to non-undulatory perturbations is retrieved by setting = in Eq, ()3.34|) 

The results of the numerical calculations that were performed for this special case are presented in Figs. [6] and [7] 
show good agreement with the analytical estimate for l) — > and k ~ 0. The general case for k^O does not allow any 
tractable analytical expressions for the eigenvalues. We therefore have to resort to numerical methods. 

To find the real part of eigenvalues for the stationary solution A, we diagonalize the matrix L numerically using 
the QR algorithm from Ref. [22] . We verified our results with the LAPACK linear algebra package |25] , and verified 
convergence of the eigenvalues up to n = 1000. 

Figure [H shows the real part of the least stable eigenvalue, i.e. that with smallest real part, as a function of k. 
The dependence on k is quadratic and for fc = the smallest eigenvalue is equal to 0. The stability decreases with 
decreasing D; for D decreasing, at a fixed value of fc = 10, from D = 10^^ to D = 10^^, Re(A) changes from 1.0 to 
about 10-3. 
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Fig. 8. The eigenvalue with smallest real part as a function of k for solution A using different values of the diffusion constant 
D. The eigenvalue with smallest real part is observed to increase quadratically with k. Moreover, the eigenvalues decrease with 
decreasing D, eventually approaching zero, a signal of neutral stability. 
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smallest 100 eigenvalues 

Fig. 9. The distribution of the first 100 eigenvalues with smallest real part for fixed fc = 3 and different values of the diffusion 
constant D. The stair case is clearly recognizable and a clear accumulation of eigenvalues occurs at Re(A) = 1 and Re(A) = 
Ac = 5. 



To investigate the distribution of eigenvalues, we plot the 100 eigenvalues with the smallest real parts for fixed 
fc = 3 in Fig. [HI One can still clearly distinguish the plateaus at Re(A) = 1 and Re(A) = Ac = 5. The plateaus become 
smaller for increasing D, as the eigenvalues increase more rapidly for larger D. 

To verify our calculations, wc have performed the same calculations for the smooth model originally introduced in 
Eq. ()2.3|) (see Appendix [Cl for details). We find that the smooth model, indeed quahtatively reproduces the numerical 
results of the piecewise toy model. 

Moreover, for this model the shear banding state with two bands is also linearly stable with respect to undulations, 
as witnessed by the positivity of the real parts of the eigenvalues in Fig. [TOl Although we report here only the results 
Ac ~ 5; we verified linear stability of the two-band solution for a large number of values of G [1, 10]. This agrees 
with the picture that arose from our analytical study for the case of zero diffusion, where wc found that one eigenvalue 
is equal to zero whereas the others are concentrated at Ac and 1. The effect of diffusion is to shorten the plateaus 
where eigenvalues have real part 1 or Ac, and this behavior persists for the entire range of Ac studied. 

The continuous (Fig. fTOl a)) and piecewise (Fig. [8]) models display qualitatively similar features: the curves initially 
increase roughly quadratically with fc and are always positive. For larger k values we find Re(A) ~ fc, which suggests 
that for this fc regime the imaginary components of the evolution of the amplitudes dominate; and for still larger values 
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k 400 smallest eigenvalues 



Fig. 10. Real part of the eigenvalue with the smallest real part, as a function of the undulation wave number k for four different 
values of the diffusion constant D for the continuous model (a). For small fc < 1 the dependence of Re(A) is again approximately 
quadratic. For larger k bumps occur, which are probably caused by the fact that the imaginary part of the equations start 
to dominate the eigenvalues. The distribution of 400 eigenvalues with smallest real part (b) for fc = 3 is very similar to the 
distribution obtained by the piecewise toy model of Fig. |9] 



of k the dependence of Re(A)(fc) levels off. This feature was occurs in the piecewise model at much larger fc-values 
not shown in Fig. [H In both cases smaller D gives rise to less stable eigenvalues. Precise comparison between the two 
models is difficult as the wavenumber k is scaled with a factor a = 20 in the continuous case. This parameter a is 
necessary in order to find a stable shear banding solution in the continuous case [3]. One striking difference between 
the curves in Fig. [8] and Fig. [TOT a) is the behavior in the limit k 0. In the smooth model with fc ^ the smallest 
eigenvalue depends on D. The eigenvalue results only when Z? ^ 0, in contrast to the piecewise model which is 
always neutrally stable with respect to non-undulatory perturbations. This can mathematically be understood from 

the fact that the interface width scales as ^J~B for D ^ Q. Therefore the matching at the interface can be done in 
a smoother fashion when D becomes increasingly smaller. In the limit D ^ Q this will result in a solution to the 
piecewise model which approximates the solution to the smooth model. If we compare Fig. [TOTb) with Fig. [9l we 
see the same qualitative behavior. Again plateaus occur, but because of the smoothness of the model the curves are 
more regular than those obtained for the piecewise toy model. Nevertheless, Fig. [TU] suggests that the results for the 
piecewise model arc generic and also valid for general smooth models without normal stresses. 

From our findings we therefore conclude that normal stresses are generically responsible for rendering the Johnson- 
Scgalman model linearly unstable for long wavelength undulations. This confirms the results which were first conjec- 
tured in Refs. [I41I20]. 



4 Discussion 



We have studied, to a large extent analytically, a toy model for shear banding without normal stresses, but with spatial 
gradient terms (stress diffusion) . We captured the general characteristics of a nonmonotonic stress-shear rate relation 
by introducing a piecewise linear stress-shear rate curve, which makes analytical calculations tractable. The results 
obtained for the piecewise model were shown to be in qualitatively agreement from those obtained using a smooth 
function g. 

For the set of parameters chosen, we found multiple stationary solutions when D < 10~^. For D = 10~^ we obtain 
three stationary states: one being the commonly observed two-band profile and the other two having three and four 
bands. The two-band profile was shown to be linearly stable with respect to two-dimensional undulations. The linear 
stability of the two-band profile is at variance with analogous results for Johnson-Segalman (JS) model, which was 
shown by Fielding to have linear instabilities for a certain range of k vectors [TH] with fc « 1. This strongly suggests 
that normal stresses, absent in our model, are responsible for the linear instabilities arising in the JS model. Related 
recent work by Fielding and co-workers [T91I20] showed that these linear instabilities are suppressed if nonlinear effects 
are taken into account. It might be that nonlinear perturbations induce instabilities in the scalar model studied here. 
The behavior of such (weakly) nonlinear instabilities is a subject of interest and future research. 
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A Integration Constants 

The integration constants ci,C2, di, c?2 can be expressed in terms of yo, the spatial position at which 7 equals 70. From 
the four equations (|2. 15112. 16|) we find 



ci 



C2 



70 



cosh (yjj/o) 

71 — aS 



cosh 



(yi - 1)/VD 



{aE - ^c7o)cos \ ^/-gyo 



tanh 



tanh ( ,/^2/o] sin [J %yo) (^c7o - aU) 



{\f§yo) cos (y^ yo) {Aclo - aS) 



(A.la) 
(A.lb) 

(A.lc) 
(A.ld) 



From these equations, together with Eqs. 



the values of yo and y\ can in principle be 



where z = 

71—70 

determined as a function of S. To fix the value of the selected stress S* an additional equation is required. From the 
global constraint of an imposed average shear rate (7) we can find an expression for S*. We substitute Equations (jA.ip 
for the constants ci^2 and cZi_2 into Eqs. (|2.14all2.14cp . and calculate the average value of the shear rate directly by 
splitting up the integral in three parts: [0, yo]j[yoj yi]i [yi, !]• By equating the sum of the three parts to (7), we obtain 
the following expression for S in terms of ya,yi, (7), and 70,71: 



Ac' 


^o{Ac + z){yo~yi)-\/^ 


^0 tanh (yjyo) (^) + ii tanh 






1) 




Ac{yo - yi + zyi) + yoz - \/d 


tanh (y^zjo) (^) + Ac{l + z) tanh 









(A.2) 



B Amplitude evolution equations 

The time evolution of the amplitudes an(fc, i) and (57o(fc, t) can be determined by substituting (57^ (x, y, i) in Eq. (|3.12ap . 
When we project out the cosine modes and use the definition of the scaled amplitudes a„(fc,t) = a„(fc,t)(2n + 
l)7r(-l)"/4, we obtain 



dam 
di 



97 



ID 



dy 



dj:6yo{k, i) - d{adS{k, t) = [A^ + Dk^) 



97 



ID 



dy 



Sya{k,i)- \ A^ + Dk^ + D 



7r(2m+ 1)' 



2yo 



+ Dk'^aS{k,i) - 

00 



ifc(-l)™(2m+ l)7r 
Wo 

4«„(-l)" 



(2n-f 1)tt 



16%' 



pA _ rpA~\ 

J- Ijiv-r-i /y-i ' 111- 



97 



ID 



(2n + l)37r3 L -run 



dyo 



Syoik,i)F,i 



dy 



5yoik,t) 



yo 



(B.l) 
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ik 

yo 



dan 



di (2n+l) 



5] a. 



^ a„(l-(-l)")16^g . ^ 
^ (2n+l)37r3 ^° ^„ 



^-4- + (70 - -4-)^ 



(2n+ l)27r2 



^ +(70-— )Vi? 



70 



= Dk^aSS - {Ac + Dk^)S% - 8{Ac + Dk'^) ^ 



n=0 ' ^ WO 



yo7o'57o 



2L> ^ ik 
^ yo 



(B.2) 



In Eq. (jB.2p wc defined the real function iTi (fc, as the imaginary contribution to the evolution of Sjq. The coefficients 
R^^,T^, F^, L^^, in Eqs. (jBTp are defined in terms of integrals of functions on the interval [0, yo]'- 



Vo 



cos 



cos 



(2m + \)Tiy\ ( (2n + l)7r?;\ 



2yo 



V 2yo 



dy 



■dy 



2yo 



{2m+l)'Ky\ f d^^ 



2 dy 



■ dy 



Vo 



cos 



{2m + l)T:y\ f {2n + l)7ry\ aSy 



2yo 



2ya J Ac 



dy, 



(B.3a) 
(B.3b) 
(B.3c) 



and = _i, and = 1. The evolution equations for the amplitudes 6n{k,i) is derived by substituting 

the expression for 6'y^^^{x, y, i) in Eq. p.l2cp . projecting out the cosine mode,s and using the scaling relation between 
dn and (5„ , yielding 



dSm 
di 



ID 



dy 



diSyi{k, i) - dia6S{k, t) = (1 + Dk^) 



ID 



Vl 



dy 



Syi{k,i) - [l + Dk^ +D 



yi 



7r(2TO-h 1) 



bk^aS{k, i) 



OC 



T I r) p 



ik{-~l)'"\2m + 1)tt 
2(1 -yi) 

4(5„(-l)' 



ri=l 



{2n + l)n 



{2n+ l)37r3 

- ^dY^ 

dy 



hD _ rfD 

mn ni 



dj 



ID 



Syi{k,i) 



dyi I , 

/ yi 

- ik{j)d, 



2(1 -yi) 

Syi{k,i)F,f, 



ik 



dj 



ID 



dy 



5yi 



(B.4) 



where the functions T^, Fj^,L^^, are defined similarly to Eqs. (jB.Sp . with yo everywhere replaced by 1 — j/i 
and Ac set to 1. If we integrate the amplitude equation for dn{k, t) from yi to 1, we obtain the evolution equation for 
the 5ji{k, i), 

n— n— \ / n— 



(i-yi) 



f. ^„(i-(^ini6(i-yi)2 , ^ 

^ (2n+l)37r3 ^i + Z^ 



8(1 - yi)dn 



ri=0 



(2n + 1)2^2 

n— ^ ^ 



(7) + ai7(yi - 1) - (71 - aU)V b 



{i)ii'm)-^^^^^^-{i.-am 



(1 - yi)'57i ((7) - aS{l - yi) + (^ji - as'j 



(1 - yi)^7i^7 



- \ =Dk^a6S - (1-f Dp)6ji - 8(1 + Dk^) ^ 



2D 



7D 
ik 



n2{k,i), 

(B.5) 
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where i72(fc,i) is the imaginary contribution to the evolution of S^i. The equation for the amplitudes &„(fc,i) can 
be found from the evolution equation of q (Eq. I3T6P and substituting expression (j3T5|) . which leads to a lengthy 
differential equation for hn{k,t). After projecting out the sine modes over [yo7yi]j we find 



dt 2 



dy 



dy 



Syo 



- -dtadU{l - (-1)™) = /3„ 



- Dk^ - D 



(2M +l)TT-h 



+ 



2 

ikmn 



dy 

' Vl 



dy 



Syo 



yo 



+ -Dk^aSU{l - {-ID 



2(yi - yo 



dy L 2 

/ Vo 



-diPi{m) + d2P2im)] + yaSjja 



97 



ID 



dy 



-d,Vi{m) + d2V2{m)] 



+ E 

00 

-E 



(2?! + iflT^ 

4/3,1 (yi - yo) 



[-diPi{m) + d2P2im)]-J2 



1/0 a,, 



n=0 



{2n + 1)^71^ 



'diVi{m) + d2V2{m) + diyoPi{m) - d2yQP2{m)] 



-diVi{m) + d2V2{m) + diyoPi{m) - d2yo-P2(m)] + ^ 



4/3n(yi - yo)^ ADPn 



X [~diGi(m,n) + d2G2(m,n)\ + } 62 H 



n=l 



97 



ID 



l-yo + (yi-l)(-iryi 



yo I I oyi - yi 



dy 



dy 



Syo 



z \{2M + f)7r - /i 



z V(2M + l)7r-/i 



[diViim) - d2V2im,)] 



[diPiim) - d2P2im)] 



dj 



ID 



dy 



Syo 



dy / , 

^ yi 

Ci ~ C2W{m) 



Syi 



97 



ID 



dy 



Syo 



Ci 



''yi - ^t:^ (jyo 



9y 



9y 



9y 



(^yi 



yofiyi - 



dy 



dy 



yiSya 



+ aC2{yi - yo)SS ^ yo + yii 1)"') _ [diPi(m) - d2P2(™)] - \ - 



- (-f)™)^i ^Ci 



(B.6) 



where 



Pi (to) = 
yi(m) = 
Gi(to, 71) 
Qi{m,n) 
W{m) = 



z 

b 

z 

b 



. ( mniy - yo) 
sm I sm 

yo V yi - yo 

. / m7r(y - yo) 
sm I I sm 

yi - yo 



yo 



yi 



I — I sm 
D 



ya 



mrjy - yo) 
yi - yo 



sm 



:y j dy 

-y^ ydy 

7r(y - yo 
yi - yo 



sm 



yi 



ya 



,mr{y-yo)\ . f rmT{y - yo" 
sm I — : z I sm | — z — I y dy 



yi - yo 



y 



Vo 



. ( mr{y - yo) 
sm I — ^— 

yi - yo) / yi - yo 



yi - yo) 

dy. 



'-^v ] ydy 



(B.7a) 
(B.7b) 
(B.7c) 
(B.7d) 
(B.7e) 
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The functions P2{m), V2{m), G'2(to, n) are obtained if the sin yy jjVj function is replaced by cos ) The constants 

Ci and C2 are defined by 



Ci = yo{z + + V 4- ho - -7- 



c?i sin ( ^1 —ijo ) — ^2 cos 



-yo 



C2 = 70 



Ac7o - aS 



(B.8a) 
(B.8b) 



The final evolution equation we need is that for SE, which is obtained by using Eqs. (|B.2|) and (|B.5p and substituting 
the expressions for dtS^o, dtSji in the differentiated Eq. (|3.18p . We then find that the evolution equations for the 
average stress perturbation SE read: 



ad,SE - A{y, - yo) ^ -^J^^^^J^ ^im - Vo) E (2n + W + '^^^ " ^ 

Tl— 



d/32 



- (2n+l)27r2 '^"^ dt (2n+l)27r2 

-0 ^ ^ n— ^ ^ 



D{2-yo-yi) 



n=0 ' 



yo ^o(2n+ 1)2^2 



E (2n+"l)27r2 + + " ^1 ~ + + ^^'^(^1 + ^"^T^ " "^^'"^^ 



+ I ni(k,t) + — ^— 7T2(fc,t) 



2yo 



2(1 -yi) 



(B.9) 



C Continuous model 

In this appendix we present the results of stability analysis of the first model, in which the nonmonotonic stress is 
encoded in the smooth function g(^) (Eq. 12. 3p . A straightforward calculation shows that for this case the amplitudes 
a„(fc, t) arc not given by Eq. (jB.ip . but instead obey 

danjk, t) 

n=0 '^^ ri=0 



E ' — cos(n7ry) = - E -D(n^7r^ + k^) a„(fc, i) cos(ri7ry) - ag'[a{{aY') + {■y)/a~ (Tp^(?/))] E an{k,i) cos{mTy) 

n—l 

(C.l) 



+ iklaY,ar,ikJ) 1 dy&l''{y)-Y,an{k,i)v^''iy)cos{mry) 



n=0 



Projecting out the cosines, we have the following matrix equation 

a = - (P™"* • a - 2iak Q"°"* • a) , (C.2) 

where the dot denotes denotes differentiation with respect to time. The entries of matrix P'^Q"* arc given as 

Pg-"* = (1 + ne) (C.3) 

PoT* = « f\'H{al'') + {^)/a-al^{y))]cos{n7Ty)dy, {n>l) (C.4) 
Jo 

P„T* = 0, (n>l) (C.5) 

P^T ^{^ + Din^T^ + e))6nra + 2a [ g'[ai{al'') + {^)/a~al°iy))]cosimry)cosim7ry)dy, (m,n>l). (C.6) 
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For the matrix Q"^""*, wc find the entries 



f^cont 

Qcont 
On 



-\cont 
InO 



cont 
nn 



Q 



cont 
nin 



•ID\ 



y^l^iy) dy 



^0 (n even) 

2((o'p°) + -^) fl ^lD^ \sin(mTy) , / i in 

Jo <^P (y) L dy (nodd) 
{n even) 



.-Io^l''(yr-^dy 



(7) 



3 
Ann 



7; [ ^l^{y)ydy-— [ sm{mry)al^ (y) dy 
^ Jo Jo 



1 - (-1)" 



~1D\ 



2 



-"(1) 



-2m2 - 2n2 + (-l)™-"(77;, + + (-1)™+"(to - n) 
27r^(m — n)^{m + n)^ 



, (">1) 



^ {y) [msin((n — m)TTy) + nsin((n — m)iTy) — msin((n -I- m)7Ty) + nsin((n + m)7ry)] 



27r^(n^ — m^) 



smijiTTy) cos(m7r?/) 



cos(n7ry) — 1 



a^-(2/)dy+(-l)"a^-(l) 



TOTT sin(m7ry) dy, (m^n, m,n>\) 



(C.7) 
(C.8) 

(C.9) 



(C.IO) 



(C.ll) 



All the matrix elements have to be computed numerically, after which the linear stability of the solutions can be 
examined along the same lines as for the piecewise model. From these computations, we conclude that the continuous 
equivalent of the piecewise system exhibits the same stability behavior as the piecewise model; it is linearly stable. For 
long wavelength undulations the dispersion relation is again dominated by the real diagonal matrix elements. That is, 
we find again a quadratic dependence of the eigenvalue with the smallest real part on k. 
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